home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / PIBSIGS / CDBETA.PAS < prev    next >
Pascal/Delphi Source File  |  1986-06-22  |  8KB  |  216 lines
  1. (*-------------------------------------------------------------------------*)
  2. (*               CDBeta  -- Cumulative Beta Distribution                   *)
  3. (*-------------------------------------------------------------------------*)
  4.  
  5. FUNCTION CDBeta(     X,     Alpha, Beta: REAL;
  6.                      Dprec, MaxIter    : INTEGER;
  7.                  VAR Cprec             : REAL;
  8.                  VAR Iter              : INTEGER;
  9.                  VAR Ifault            : INTEGER  ) : REAL;
  10.  
  11. (*-------------------------------------------------------------------------*)
  12. (*                                                                         *)
  13. (*       Function:  CDBeta                                                 *)
  14. (*                                                                         *)
  15. (*       Purpose:   Evaluates CPDF of Incomplete Beta Function             *)
  16. (*                                                                         *)
  17. (*       Calling Sequence:                                                 *)
  18. (*                                                                         *)
  19. (*            P     := CDBeta(     X, Alpha, Beta: REAL;                   *)
  20. (*                                 Dprec, Maxitr : INTEGER;                *)
  21. (*                             VAR Cprec         : REAL;                   *)
  22. (*                             VAR Iter          : INTEGER;                *)
  23. (*                             VAR Ifault        : INTEGER  ) : REAL;      *)
  24. (*                                                                         *)
  25. (*                 X      --- Upper percentage point of PDF                *)
  26. (*                 Alpha  --- First shape parameter                        *)
  27. (*                 Beta   --- Second shape parameter                       *)
  28. (*                 Dprec  --- Number of digits of precision required       *)
  29. (*                 Maxitr --- Maximum number of iterations                 *)
  30. (*                 Cprec  --- Actual resulting precision                   *)
  31. (*                 Iter   --- Iterations actually used                     *)
  32. (*                 Ifault --- error indicator                              *)
  33. (*                            = 0:  no error                               *)
  34. (*                            = 1:  argument error                         *)
  35. (*                                                                         *)
  36. (*                 P      --- Resultant probability                        *)
  37. (*                                                                         *)
  38. (*       Calls:                                                            *)
  39. (*                                                                         *)
  40. (*            ALGama                                                       *)
  41. (*                                                                         *)
  42. (*       Method:                                                           *)
  43. (*                                                                         *)
  44. (*            The continued fraction expansion as given by                 *)
  45. (*            Abramowitz and Stegun (1964) is used.  This                  *)
  46. (*            method works well unless the minimum of (Alpha, Beta)        *)
  47. (*            exceeds about 70000.                                         *)
  48. (*                                                                         *)
  49. (*            An error in the input arguments results in a returned        *)
  50. (*            probability of -1.                                           *)
  51. (*                                                                         *)
  52. (*-------------------------------------------------------------------------*)
  53.  
  54. VAR
  55.    Epsz : REAL;
  56.    A    : REAL;
  57.    B    : REAL;
  58.    C    : REAL;
  59.    F    : REAL;
  60.    Fx   : REAL;
  61.    Apb  : REAL;
  62.    Zm   : REAL;
  63.    Alo  : REAL;
  64.    Ahi  : REAL;
  65.    Blo  : REAL;
  66.    Bhi  : REAL;
  67.    Bod  : REAL;
  68.    Bev  : REAL;
  69.    Zm1  : REAL;
  70.    D1   : REAL;
  71.    Aev  : REAL;
  72.    Aod  : REAL;
  73.  
  74.    Ntries : INTEGER;
  75.  
  76.    Qswap  : BOOLEAN;
  77.    Qdoit  : BOOLEAN;
  78.    Qconv  : BOOLEAN;
  79.  
  80. LABEL   20, 9000;
  81.  
  82. BEGIN (* CdBeta *)
  83.  
  84.                                    (* Initialize *)
  85.    IF Dprec > MaxPrec THEN
  86.       Dprec := MaxPrec
  87.    ELSE IF Dprec <= 0 THEN
  88.       Dprec := 1;
  89.  
  90.    Cprec  := Dprec;
  91.  
  92.    Epsz   := PowTen( -Dprec );
  93.  
  94.    X      := X;
  95.    A      := Alpha;
  96.    B      := Beta;
  97.    QSwap  := FALSE;
  98.    CDBeta := -1.0;
  99.    Qdoit  := TRUE;
  100.                                    (* Check arguments *)
  101.                                    (* Error if:       *)
  102.                                    (*    X <= 0       *)
  103.                                    (*    A <= 0       *)
  104.                                    (*    B <= 0       *)
  105.  
  106.    Ifault := 1;
  107.  
  108.    IF( X <= 0.0 ) THEN GOTO 9000;
  109.  
  110.    IF( ( A <= 0.0 ) OR ( B <= 0.0 ) ) THEN GOTO 9000;
  111.  
  112.    CDBeta := 1.0;
  113.    Ifault := 0;
  114.                                    (* If X >= 1, return 1.0 as prob *)
  115.  
  116.    IF( X >= 1.0 ) THEN GOTO 9000;
  117.  
  118.                                    (* If X > A / ( A + B ) then swap *)
  119.                                    (* A, B for more efficient eval.  *)
  120.  
  121.    IF( X > ( A / ( A + B ) ) ) THEN
  122.       BEGIN
  123.          X      := 1.0 - X;
  124.          A      := Beta;
  125.          B      := Alpha;
  126.          QSwap  := TRUE;
  127.       END;
  128.  
  129.                                    (* Check for extreme values *)
  130.  
  131.    IF( ( X = A ) OR ( X = B ) ) THEN GOTO 20;
  132.    IF( A = ( ( B * X ) / ( 1.0 - X ) ) ) THEN GOTO 20;
  133.    IF( ABS( A - ( X * ( A + B ) ) ) <= Epsz ) THEN GOTO 20;
  134.  
  135.          C     := ALGama( A + B ) + A * LN( X ) +
  136.                   B * LN( 1.0 - X ) - ALGama( A ) - ALGama( B ) -
  137.                   LN( A - X * ( A + B ) );
  138.  
  139.          IF( ( C < -36.0 ) AND QSwap ) THEN GOTO 9000;
  140.  
  141.             CDBeta := 0.0;
  142.             IF(  C < -180.0 ) THEN GOTO 9000;
  143.  
  144.                                    (*  Set up continued fraction expansion *)
  145.                                    (*  evaluation.                         *)
  146.  
  147. 20:
  148.       Apb    := A + B;
  149.       Zm     := 0.0;
  150.       Alo    := 0.0;
  151.       Bod    := 1.0;
  152.       Bev    := 1.0;
  153.       Bhi    := 1.0;
  154.       Blo    := 1.0;
  155.  
  156.       Ahi    := EXP( ALGama( Apb ) + A * LN( X ) +
  157.                      B * LN( 1.0 - X ) - ALGama( A + 1.0 ) -
  158.                      ALGama( B ) );
  159.  
  160.       F      := Ahi;
  161.  
  162.       Iter   := 0;
  163.  
  164.                                    (* Continued fraction loop begins here. *)
  165.                                    (* Evaluation continues until maximum   *)
  166.                                    (* iterations are exceeded, or          *)
  167.                                    (* convergence achieved.                *)
  168.  
  169.       Qconv  := FALSE;
  170.  
  171.       REPEAT
  172.  
  173.          Fx     := F;
  174.  
  175.          Zm1    := Zm;
  176.          Zm     := Zm + 1.0;
  177.          D1     := A + Zm + Zm1;
  178.          Aev    := -( A + Zm1 ) * ( Apb + Zm1 ) * X / D1 / ( D1 - 1.0 );
  179.          Aod    := Zm * ( B - Zm ) * X / D1 / ( D1 + 1.0 );
  180.          Alo    := Bev * Ahi + Aev * Alo;
  181.          Blo    := Bev * Bhi + Aev * Blo;
  182.          Ahi    := Bod * Alo + Aod * Ahi;
  183.          Bhi    := Bod * Blo + Aod * Bhi;
  184.  
  185.          IF ABS( Bhi ) < Rsmall THEN Bhi := 0.0;
  186.  
  187.          IF( Bhi <> 0.0 ) THEN
  188.             BEGIN
  189.                F     := Ahi / Bhi;
  190.                Qconv := ( ABS( ( F - Fx ) / F ) < Epsz );
  191.             END;
  192.  
  193.          Iter  := Iter + 1;
  194.  
  195.       UNTIL ( ( Iter > MaxIter ) OR Qconv ) ;
  196.  
  197.                                    (* Arrive here when convergence    *)
  198.                                    (* achieved, or maximum iterations *)
  199.                                    (* exceeded.                       *)
  200.  
  201.    IF ( Qswap ) THEN
  202.       CDBeta := 1.0 - F
  203.    ELSE
  204.       CDBeta := F;
  205.  
  206.                                    (* Calculate precision of result *)
  207.  
  208.    IF ABS( F - Fx ) <> 0.0 THEN
  209.       Cprec := -LogTen( ABS( F - Fx ) )
  210.    ELSE
  211.       Cprec := MaxPrec;
  212.  
  213. 9000:  (* Error exit *)
  214.  
  215. END   (* CDBeta *);
  216.